# SI Figure 6
rm(list=ls())
gc()

quartz(width = 10, height = 5, type = "pdf", "SI Figure joseph", file = "SI_Figure_joseph.pdf")
par (fig=c(0,1,0,1), # Figure region in the device display region (x1,x2,y1,y2)
	omi=c(.8,0,0,.5), # global margins in inches (bottom, left, top, right)
	mai=c(0,.8,1,.3)) # subplot margins in inches (bottom, left, top, right)
	layout(matrix(c(1,2), 1, 2, byrow = TRUE), widths=c(3,3), heights=c(3,3))

# Panel 1:
t1 <- data.frame(r1=c(rep(1,7),rep(2,7)), r2=c(rep(2,7),rep(1,7)))
plot(1:10,1:10,type="n",xlim=c(0,14.5),ylim=c(0,3),xaxs="i",yaxs="i",bty="L",yaxt="n",xaxt="n",xlab="",ylab="")
lines(1:14, t1$r1)
points(1:14, t1$r1, pch=16)
lines(1:14, t1$r2)
points(1:14, t1$r2, pch=16)

# Panel 2:
t2 <- data.frame(r1=c(1,1,2,1,1,2,2,2,1,2,2,1,1,2),r2=c(2,2,1,2,2,1,1,1,2,1,1,2,2,1))
plot(1:10,1:10,type="n",xlim=c(0,14.5),ylim=c(0,3),xaxs="i",yaxs="i",bty="L",yaxt="n",xaxt="n",xlab="",ylab="")
lines(1:14, t2$r1)
points(1:14, t2$r1, pch=16)
lines(1:14, t2$r2)
points(1:14, t2$r2, pch=16)


mtext(side=1, expression("Time"), outer=T, line=1, cex=.9, at=0.4)
mtext(side=1, expression("Time"), outer=T, line=1, cex=.9, at=0.9)
mtext(side=3, expression(Delta * "T"), outer=T, line=-4, cex=.9, at=0.1)
mtext(side=3, expression(Delta * "T"), outer=T, line=-4, cex=.9, at=0.6)


dev.off()




# Temperatures: Equal means and variances
sum(t1$r1) + sum(t1$r2) == sum(t2$r1) + sum(t2$r2) # Average temperatures
sum(t1$r1) - sum(t1$r2) == sum(t2$r1) - sum(t2$r2) # Variance temperatures

# Damages/Consumption: Equal means and variances
d <- function(x) { 1/(1+.0023*x^2) }
sum(d(t1$r1)) + sum(d(t1$r2)) == sum(d(t2$r1)) + sum(d(t2$r2)) # Average damages
sum(d(t1$r1)) - sum(d(t1$r2)) == sum(d(t2$r1)) - sum(d(t2$r2)) # Variance damages

# Utility: Equal means and variances
u <- function(x) { x^(-.5)/(-.5) }
sum(u(d(t1$r1))) + sum(u(d(t1$r2))) == sum(u(d(t2$r1))) + sum(u(d(t2$r2))) # Average utility
sum(u(d(t1$r1))) - sum(u(d(t1$r2))) == sum(u(d(t2$r1))) - sum(u(d(t2$r2))) # Variance utility

# With growth: Equal means but different variances
g <- 1.1^(0:13)
sum(d(t1$r1)*g) + sum(d(t1$r2)*g) == sum(d(t2$r1)*g) + sum(d(t2$r2)*g) # Average damages
sum(d(t1$r1)*g) - sum(d(t1$r2)*g) == sum(d(t2$r1)*g) - sum(d(t2$r2)*g) # Variance damages

sum(u(d(t1$r1))*g) + sum(u(d(t1$r2))*g) == sum(u(d(t2$r1))*g) + sum(u(d(t2$r2))*g) # Average damages
sum(u(d(t1$r1))*g) - sum(u(d(t1$r2))*g) == sum(u(d(t2$r1))*g) - sum(u(d(t2$r2))*g) # Variance damages

# With discounting
df <- 1/(1.1^(0:13))
sum(d(t1$r1)*df) + sum(d(t1$r2)*df) == sum(d(t2$r1)*df) + sum(d(t2$r2)*df) # Average damages
sum(d(t1$r1)*df) - sum(d(t1$r2)*df) == sum(d(t2$r1)*df) - sum(d(t2$r2)*df) # Variance damages

# With growth and discounting
sum(d(t1$r1)*g*df) + sum(d(t1$r2)*g*df) == sum(d(t2$r1)*g*df) + sum(d(t2$r2)*g*df) # Average damages
sum(d(t1$r1)*g*df) - sum(d(t1$r2)*g*df) == sum(d(t2$r1)*g*df) - sum(d(t2$r2)*g*df) # Variance damages

# Any time-dependent transformation (growth, discounting) affects the variance of NPV damages, but not the mean.
# It will also increase the variance of u, but not the mean, and therefore not the risk premium.
